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The recently proposed full configuration interaction quantum Monte Carlo method allows access to essentially 
exact ground-state energies of systems of interacting fermions substantially larger than previously tractable 
without knowledge of the nodal structure of the ground-state wave function. We investigate the nature of the 
sign problem in this method and how its severity depends on the system studied. We explain how cancelation 
of the positive and negative particles sampling the wave function ensures convergence to a stochastic repre- 
sentation of the many-fermion ground state and accounts for the characteristic population dynamics observed 
in simulations. 



I. INTRODUCTION 

One of the major goals of electronic structure methods 
is to produce accurate ground-state energies and prop- 
erties of many-electron systems.^ Quantum chemistry 
provides a hierarchy of ah initio methods^ based upon 
Hartree-Fock, of which coupled-cluster singles and dou- 
bles with perturbative triples (CCSD(T)) is the most ac- 
curate applicable to medium-sized molecules.^ Quantum 
Monte Carlo methods such as Green's function Monte 
Carlo^i^ (GFMC), the closely related diffusion Monte 
Carlcfi"— (DMC), and auxiliary- field quantum Monte 
Carlo^ (AFQMC) produce accurate results via a stochas- 
tic sampling of the many-electron wave function, but 
none of these methods is exact: GFMC and DMC sim- 
ulations of all but the smallest systems converge to the 
physically irrelevant many-boson ground state unless the 
fixed- node approximation is made^ii whilst AFQMC re- 
quires the phaseless approximation^^ in order to avoid 
an exponential growth in noise except in certain special 
casesJi 

The full configuration interaction (FCI) method^^ 
casts the Schrodinger equation as a matrix eigenvalue 
problem, in which the requirement that the many- 
electron wave function be anti-symmetric with respect to 
exchange of electrons is imposed by working in a space 
of Slater determinants formed from a finite basis set of 
single-particle wave functions. The lowest eigenvalue and 
eigenvector of the FCI Hamiltonian matrix give the ex- 
act ground-state energy and wave function of the sys- 
tem, subject only to the error due to the finite basis set. 
Whilst the computational cost of FCI scales factorially 
with system size, it nevertheless represents the holy grail 
of electronic structure methods. 

In 2009, Booth, Thom and Alavi^^ introduced a new 
stochastic approach in which the nodal structure of the 
ground-state wave function emerges spontaneously by 
samphng the discrete space of Slater determinants. Their 
"full configuration interaction quantum Monte Carlo" 
(FCIQMC) method yields exact (i.e. FCI-quality) results 
whilst requiring a fraction of the memory of an FCI cal- 



culation using the same basis. The memory required by 
FCIQMC simulations still scales factorially with system 
size, but the exponent appears to be substantially smaller 
than for FCI simulations. Moreover, unlike the iterative 
diagonaliation schemes required for FCI, the FCIQMC 
algorithm is readily parallelizable and can run efficiently 
on thousands of cores. Alavi and co-workers have used 
FCIQMC to reproduce essentially every molecular FCI 
calculation ever done and have obtained ground-state en- 
ergies for systems with Hilbert spaces many orders of 
magnitude larger than the largest FCI calculations^^ — 
We believe that the FCIQMC method will become in- 
creasingly important in the electronic structure commu- 
nity, especially if improved algorithms or substantially 
cheaper approximations can be found. Our motivation 
for exploring the behavior of the method is to provide 
insight into possible improvements. 

An FCI calculation based on iterative diagonalization 
(using, for example, the Davidson method) requires the 
storage of at least two vectors, each of length equal to the 
size of the Hilbert space. An FCIQMC simulation using 
the same basis requires the storage of the labels of the de- 
terminants occupied by a population of stochastic walk- 
ers (which, following Anderson,^ we call "psi-particles" 
or psips) scattered over the same Hilbert space. In order 
for FCIQMC to be more efficient than FCI, the num- 
ber of psips required must be a small fraction of the 
size of the Hilbert space. The fraction required is sys- 
tem dependent and provides a measure of how "hard" it 
is for FCIQMC to treat that system. The hardness is 
surprisingly difficult to predict. FCIQMC is wildly suc- 
cessful for some systems, such as the neon atom, where 
it requires only ~0.01% of the memory of a conven- 
tional FCI calculation Even for the nitrogen molecule, 
a classic example of a strongly correlated system and 
a tough test for quantum chemical methods, FCIQMC 
used only a quarter of the memory of the equivalent FCI 
calculation.^^ However, FCIQMC struggles to describe 
the methane molecule for which Hartree-Fock is a very 
good approximation. We show in Sec. [Ill that FCIQMC 
also struggles when applied to the Hubbard model and 
cannot treat systems larger than existing FCI methods 
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unless U is very small. What is it that makes a system 
difficult? Evidently the answer is more complicated than 
whether or not the system is strongly correlated. 

The aim of this paper is to understand the FCIQMC 
algorithm better. Why are some systems more difficult 
to treat than others? What determines the characteristic 
population dynamics observed in FCIQMC simulations? 
How does the cancelation of positive and negative psips 
ensure convergence to the many-fermion ground state? 
How many psips are required to obtain correct results? 
What goes wrong if the population of psips is too small? 
We provide at least partial answers to all of these ques- 
tions. 

Section |TT1 reviews the FCIQMC method and applies 
it to the Hubbard model as an example. In Sec. ITTTl 
we discuss the nature of the sign problem in FCIQMC 
and explain the effect of canceling positive and negative 
psips. Section llVl shows how our analysis of the sign prob- 
lem also explains the characteristic population dynamics 
observed in FCIQMC simulations. Section |Vl describes 
several special cases in which the sign problem can be 
removed entirely. We offer some concluding remarks in 
Sec. ED 



and impose the normalization constraint using a La- 
grange multiplier, the optimal coefficients coi form the 
lowest eigenvector cq of the matrix-eigenvalue problem 



(5) 



where Eq is both the Lagrange multiplier and the ground- 
state energy eigenvalue. This approach yields the FCI 
wave function for the finite basis set used. 

One way to find the ground-state eigenvector cq of 
the Hamiltonian matrix H is to solve the imaginary-time 
Schrodinger equation: 



(6) 



The solution vector c(r) converges to cq as r ^ oo when- 
ever the starting vector c(r=0) has a non-zero overlap 
with Co; indeed, this is also the driving principle behind 
other projector methods such as DMC and GFMC. The 
convergence is easy to demonstrate by considering the 
formal solution of the imaginary-time Schrodinger equa- 
tion: 



II. FCIQMC METHOD 

We briefly summarize the FCIQMC method here, 
largely following the derivations given by Booth et al.^ 
with attention paid to details relevant later in this paper. 

Consider an orthonormal set of 2M single-particle 
spin-orbit als, {(/>!, ^2, • • • , 02m}- Previous FCIQMC 
worki^"— has used a basis of Hartree-Fock spin-orbitals, 
which is often a sensible choice but not required. One 
can construct an TV-electron Slater determinant, by 
selecting any N spin-orbitals (assuming N < 2M): 



(1) 



(2) 



To ensure that all such determinants are unique (up to a 
sign), it is convenient to require that ii < 12 < ■ - ■ < iN - 
The ground-state wave function, may be defined 
variationally as the A^-electron wave function that mini- 
mizes the energy expectation value. 



(3) 



subject to the normalization constraint = 1. If we 

restrict the search for \l/o to the subset of wave functions 
that can be expanded in the basis oi'^^CN determinants, 



(4) 



c(t) 



-Hr 



c(0). 



(7) 



If the initial vector, c(0), is expanded in terms of the com- 
plete orthonormal set of the '^^Cn eigenvectors, {cq,}, of 
the Hamiltonian matrix, H, 



c(0) =^^c.(0)Cc., 



(8) 



the solution of the imaginary-time Schrodinger equation 
can be written as 

c(r) = e-^"^^«(0)c« = ^^c.(0)e-^-"c,. (9) 



The summation is dominated by the ground-state contri- 
bution in the limit r ^ oo and thus 



c(r oo) 



^o(0)e-^°"co. 



(10) 



The steady change in normalization due to the e~^^^ 
factor is awkward but can be removed by choosing the 
zero of energy such that Eq = {). In practice, Eq is not 
known until the end of the simulation, so we instead solve 



(11) 



where the energy shift, is adjusted slowly to keep the 
normalization more or less constant. The energy shift 
converges to the ground-state eigenvalue, £^o, in the long- 
time limit. Note that previous worl^ — also subtracted 
the Hartree-Fock energy, ^hf, from the diagonal ele- 
ments of the Hamiltonian, but this amounts merely to 
a redefinition of 6*. To make the subsequent notation as 
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simple as possible, we define a "transition matrix", T, 



T = -(H-5I). 



(12) 



Note that the minimal eigenvalue of H corresponds to 
the maximal eigenvalue of T. 

The FCIQMC algorithm proposed by Booth et alM 
may be summarized as follows. Consider a collection of 
markers distributed over the space of determinants, {D\]. 
The markers do not move through the Hilbert space, so 
following Anderson^ we call them "psi particles" or psips 
instead of "walkers". (It is easy to devise alternative 
FCIQMC algorithms in which the psips do move and be- 
have much like walkers in DMC^ but we use the approach 
of Booth, Thom and Alavi here.) Each psip has both a 
location, i, and a "charge", q = ±1, associated with it. 

In one time step, Ar, we loop over the population of 
psips and allow each to "spawn" children (new psips) 
located on new determinants (which may be the same 
as the parent's determinant) according to the following 
rules: 

• The probability that a psip at i successfully spawns 
a child at j is |Tji|Ar. 

• If the spawning event is successful, the child psip 
has charge ^chiid = sign(Tji)(7parent, where (^parent is 
the charge of the parent. 

At the end of the time step, pairs of psips of opposite 
charge on the same determinant cancel each other out 
( "annihilate" ) and are removed from the simulation. For 
example, psips on determinant D\ with negative diago- 
nal elements T\\ may generate children which immedi- 
ately annihilate with the parent psips. Thus, although 
psips on different determinants may have different signs, 
all psips on any given determinant have the same sign 
at the end of every time step. Any psips which remain 
(whether originally "child" or "parent") are merged and 
are allowed to spawn new psips in subsequent time steps. 

The FCIQMC algorithn>i^ actually yields a stochastic 
sampling of the solution of the iterative equation 



Ci(r + AT) = ^(5ij+IijAT)cj(r), 



(13) 



which may be regarded as a finite- difference approxima- 
tion to the imaginary-time Schrodinger equation 



(14) 



Assuming that S = as is the case once the sim- 
ulation has converged on the ground state, the domi- 
nant eigenvectors of 5ij + ^ijAr and Tij are identical if 
Ar(£^max — £^o) < 2, where £^max is the largest eigenvalue 
of the FCI Hamiltonian matrix H. This means that there 
is no time-step error if Ar is small enough. In order for 
there to be a finite time step, the Hamiltonian must be 



bounded from both above and below. It also follows that 
the time step is limited if high-energy states are included 
in the basis; for example increasing the basis set in calcu- 
lations of quantum chemical systems requires a decrease 
in the time step. 

Rather than considering all possible spawning events, 
it is computationally efficient to allow a psip at D{ to at- 
tempt to spawn only onto its own determinant and a sin- 
gle connected determinant per time step. (Connected 
determinants are linked by non-zero transition matrix 
elements Tji.) The connected determinant is chosen 
stochastically according to some generation probability 
Pgen(j ^ i), which must be non-zero whenever Tji is non- 
zero and normalized such that XljPgen(j ^ i) = 1. To 
maximise the probability of successful spawning, the gen- 
eration probabilities should be as similar to |Tji| 
as possible; in practice, we go for speed and set the gen- 
eration probabilities for spawning on all connected de- 
terminants equal. Once the candidate spawning location 
j has been chosen, the spawning event is accepted with 
probability |Tji| Ar/pgen(j ^ i)J^ The effective size of 
the Hilbert space can be reduced by only selecting can- 
didate spawning locations of some desired symmetry. 

The population evolution of the psips over the course 
of an FCIQMC simulation has some commonly observed 
features, as illustrated in Fig.[TJ we develop simple math- 
ematical models to understand these results in Sections 
Uni and HVl The energy shift S is initially set to a value 
greater than (i.e., less negative than) the ground-state 
energy eigenvalue, ensuring that the largest eigenvalue 
of T = S*! — H is positive and hence that total psip 
population grows exponentially. Once the population 
reaches a critical size the simulation spontaneously en- 
ters a "plateau" regime, where the total psip population 
remains stable and during which the ground-state wave 
function emerges. The height of the plateau relative to 
the size of the space of Slater determinants provides a 
measure for how hard a system is to solve using FCIQMC. 
After a system- and ^'-dependent waiting time, the sim- 
ulation exits the plateau phase and the psip population 
starts to grow in an exponential fashion again, albeit at 
a slower rate than before. At this point the distribution 
of psips is a stochastic representation of the FCI ground- 
state wave function. Allowing the psip population to in- 
crease further serves only to reduce the statistical noise, 
so it is convenient to hold the population constant by 
varying the shift from this point in the simulation. The 
changes in shift must be carried out smoothly and slowly 
to avoid introducing a bias. Following Booth, we use 
the shift-update algorithm originally proposed for DMC 
simulations by Umrigar—: 



S{r + AAr) = S{r) 



AAr^n N,{r) ' ' ^^^^ 



where A is the number of iterations between which the 
shift is updated, is a damping parameter and Np{r) is 
the total psip population at time r. In the simulations 
presented here we set A = 20 and ^ = 0.1. Once the psip 
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population has settled down, the expected value qi of the 
sum of the charges of the psips on Di is proportional to 
the weight of Di in the exact ground-state wave function: 



may also be written in reciprocal space, 



lim <7i(r) oc coi- 



(16) H = J2eA^ 



Note that annihilation of psips of opposite charge does 
not affect qj. 

There are two main ways of obtaining an estimate of 
the ground-state energy from an FCIQMC simulation. 
The shift estimator is the mean value of the shift S in the 
constant-population-mode simulation after the plateau. 
The projected estimator is obtained by noting that, for 
any single determinant Dq with a non-zero ground-state 
component, the quantity 



E = lim 



{Do\H\<b{T)) 



qo 



(17) 
(18) 



tends to the ground-state energy as r ^ oo. Since the 
Hamiltonian operator only links determinants differing 
by two or fewer excitations, the number of terms included 
in the sum is limited. Note that it is important to average 
qj{r) and qo{r) separately; the projected estimator and 
its associated error can then be found by taking the ra- 
tio of the averages and using the covariance, respectively. 
Taking Dq to be the determinant with the largest overlap 
with the exact ground-state wave function minimises the 
relative stochastic noise in the denominator of the above 
equation. Furthermore, such a determinant will typically 
have single and double excitations which also have sig- 
nificant contributions to the ground-state wave function, 
and hence determinants contributing to the numerator 
will also often have a significant population. The de- 
terminant with greatest overlap may not necessarily be 
known a priori (or even be clearly defined, as is the case 
in systems studied in Sec.|Vj) but in practice the Hartree- 
Fock determinant is usually a good choice. The choice 
of 1^0 can be changed during the course of a simulation 
if a determinant with a particularly large population is 
found. 

We can illustrate the main features of the FCIQMC 
method by applying it to the Hubbard Hamiltonian^'' 
for a d-dimensional cubic lattice, 



H 



-t 



E 

(r,r') ,cr 



(19) 



where cj ^ {Cj. ^) creates (destroys) an electron of spin a 
on site r, the number operator n^. ^ = ^c^. ^ counts the 
spin-cr electrons on site r, and the summation over (r, r') 
includes all nearest-neigbor pairs of lattice sites. Assum- 
ing a system of M sites and 2M spin-orbitals subject to 
periodic boundary conditions, the Hubbard Hamiltonian 



k,cr 



U_ 

M 



^ki,t^k2,i^k3,i^ki+k2-k3,t' 



where c 



,/a7 Z^r 



ki,k2,k3 

(20) 

^ — —j= /_^^ ^^^^^'^ is the creation opera- 
tor for a Bloch state of wave vector k, the sums are 
over the M wave vectors in the first Brillouin zone, 
ek = — 2t X^iLi cos(/cia) is the one-electron kinetic energy 
for wave vector k = (/ci, /c2, • • • , ^d), and a is the lattice 
parameter (which is set to unity from now on) . The com- 
bination of wave vectors ki + k2 — ka is assumed to have 
been reduced into the first Brillouin zone by addition of 
a reciprocal lattice vector if necessary. 

Whilst the FCIQMC algorithm may be used in the 
real- and reciprocal-space representations, here, for illus- 
trative purposes, we focus on the reciprocal-space formu- 
lation, as would be appropriate \iU/t were small. As 
U /t increases and so electrons become increasingly more 
localised, the real-space formulation becomes more ap- 
propriate. Indeed, as we show in Sec. the real-space 
basis is substantially more favourable for non-zero values 
oiU/t in the one-dimensional Hubbard model. 

The imaginary-time evolution of the psip population 
and energy estimators during a k-space FCIQMC sim- 
ulation of the 2D 18-site half-filled Hubbard model at 
U/t = 4 (which is not small) is shown in Fig.[TJ As might 
be expected, the k-space FCIQMC method becomes less 
efficient at finding the ground state of this system as [//t 
increases — in fact, the psip population at the plateau 
increases approximately linearly with U/t for U/t > 2 
(Fig. [2j). When U/t > 4, one needs at least as many 
psips as there are determinants in the Hilbert space and 
the memory requirements are comparable to those of it- 
erative diagonalization. As a result, unless U/t is small, 
k-space FCIQMC is not able to treat half-filled Hubbard 
model systems substantially larger than those accessible 
to the FCI method. 

Working in a space of Slater determinants has two 
main advantages. As the basis is anti-symmetric with 
respect to exchange of two spin-orbitals, there is no need 
to use the fixed-node approximation to prevent collapse 
to the bosonic ground state, as in DMC and GFMC. Al- 
though FCIQMC still has a sign problem (see Sec. IIIip . 
the instability is not with respect to the a bosonic state 
and is often less severe. The second advantage is that 
the annihilation of psips with opposite charges proves 
surprisingly efficient in the discrete space of Slater de- 
terminants. Walker cancelation can in principle cure the 
sign problem in continuum DMC and GFMC simulations 
too^ — but such algorithms are substantially more com- 
plicated and less successful than the sign cancelation in 
FCIQMC. 
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FIG. 1. Population dynamics and energy estimators for the 
2D 18-site square lattice half- filled Hubbard model at [/ = 4t. 
Periodic boundary conditions are applied to a 45° tilted 
square cell with sides 3v^- The ground-state wave function 
has momentum k = (0,0) and the Hilbert space formed from 
all determinants of this momentum contains 1.3 x 10^ states. 
Graph (a) shows how the psip population evolves with simula- 
tion time. The psip population initially grows exponentially 
before reaching a plateau, during which the distribution of 
psips settles to model the FCI wave function. The psip pop- 
ulation then begins to grow exponentially again, reinforcing 
the wave function signal, at which point the shift is allowed 
to vary to stabilize the number of psips. Graph (b) shows the 
shift and projected energy estimators as a function of simula- 
tion time. The FCI energy is from Ref . [Tsl. FCIQMC requires 
roughly the same number of psips as the size of the Hilbert 
space to converge on the ground-state wave function in this 
case. 



coupled differential equations: 
dnt 



(21) 



These can be decoupled by adding and subtracting: 



dr 



d{nf 



dr 



(22) 



As r oo, n+ + n~ tends to the eigenvector core- 
sponding to the largest eigenvalue of T+ + T~, whilst 
n+ — n~ tends to the eigenvector corresponding to the 
largest eigenvalue of T+ — T~. We wish to find the 
latter state. However, as explained below, the largest 
eigenvalue of T+ + T~ is always larger than that of 
T+ — T~. Thus, the signal n+ — n~ always decays relative 
to n+ + n~ . Fig. [3] shows that performing an FCIQMC 
simulation without annihilation does indeed give this un- 
desired state. We note that a similar analysis has been 
performed previously for the world-line Quantum Monte 
Carlo method^ii^. 

One can prove that the largest eigenvalue of T+ + T~ 
is always greater than or equal to the largest eigenvalue 
of T+ — T~ as follows. Suppose that cq is the eigenvector 
corresponding to the largest eigenvalue, AJJJ^^, of T+ — 
T~. If Co is normalized and the Hamiltonian matrix is 
real, all components of cq may also be chosen real and so 



(23) 



III. THE ORIGIN OF THE SIGN PROBLEM IN FCIQMC 



Booth et al}^ showed that annihilation is crucial in 
FCIQMC; without it the simulation never converges to 
the FCI ground state. In this section we attempt to ex- 
plain why annihilation is required and how it helps the 
ground state to emerge. 

First, consider the effect of removing annihilation from 
the simulation procedure. This means that psips of both 
charges are permitted to reside on and propagate from 
the same determinant at the same time. We shall use 
nf and nr to represent the populations of positive and 
negative psips on determinant D\. It is convenient to 
write the transition matrix T as T+ — T~, where T+ 
contains the positive transition matrix elements, = 
max(Tjj,0), and T~ contains the absolute values of the 
negative elements, Tj7 = max(— 7^j,0). All elements of 
and T~ are therefore non-negative. The populations 
of positive and negative psips evolve according to the 



Now consider the vector |co| with components |coi|. This 
vector is also normalized and may be used as a trial state 
in the (upside-down) variational principle for the largest 
eigenvalue, A^^, of T+ + T": 



A max 
^sum 



(24) 



The right-hand side of Eq. [231 is manifestly greater than 
or equal to the right-hand side of Eq.[23l so > 

The discrete annihilation process is difficult to model 
exactl}^ in a differential equation. We therefore consider 
a simpler annihilation process in which pairs of psips of 
opposite sign on the same determinant annihilate each 
other at a constant rate 2i<i^ where is a small positive 
constant and the factor of 2 is introduced solely for alge- 
braic convenience. This leads to the differential equations 

dnf 
~d^ 



dn.^ 
~d^ 



E 

j 



I) 



(25) 
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FIG. 2. Graphs (a)-(d) show the psip population dynamics for the 2D 18-site square lattice half-filled Hubbard model at 
k = (0,0) for various values of U/t. Periodic boundary conditions are applied to a 45° tilted square cell with sides 3a/2- The 
data for U = At in (b) is identical to that in Fig. [1] Graph (e) shows that the number of psips at the plateau increases linearly 
with U/t when U /t>2. The standard errors in the numbers of psips during the plateau phases were obtained using a blocking 
analysis^^ and the straight line was fitted using the method of least squares as implemented in Ref. [l^. Unless shown, the error 
bars are smaller than the markers. 




Eq. [22] thus become 



100 150 
r/a.u. 



250 



FIG. 3. The impact of removing annihilation from the 
FCIQMC algorithm applied to a 100 x 100 randomly- 
generated real symmetric Hamiltonian matrix H. is the 
exact ground-state eigenvalue and E{t) the projected estima- 
tor for H. Ao and A(t) are the analogous quantities for the 
matrix + H~, where H = — H~ and all elements of 
and H~ are negative (i.e., all elements of and T~ are 
positive) . The simulation converges to the lowest eigenvector 
of + H~, which does not correspond to an eigenvector of 
the Hamiltonian matrix. 



d{r 



dr 



dinf 



dr 



E(^T 

j 



E 

j 



^7 



^7 



(< 
(< 



Anntn. 



The decoupled ordinary differential equations (ODEs) in 



(26) 

It is clear how annihilation enables the FCIQMC method 
to converge upon the ground state of the Hamiltonian 
matrix: as the total psip population, ^-^{nf +nj~), rises, 
the rate of annihilation events rises quadratically. This 
causes the rate of growth of the n++n~ signal to decrease 
until spawning of new psips is balanced by annihilation. 
The growth of n+ — n~, the desired solution, is not af- 
fected and so this signal eventually emerges. As shown 
in Sec. lIVI the rates of growth of the two signals actually 
become the same but due to annihilation all psips on the 
same determinant have the same sign. 

To summarize, the sign problem in FCIQMC originates 
from the in-phase combination of positive and negative 
psips, which grows at a rate determined by the largest 
eigenvalue of T+ + T~ . This eigenvalue is greater than 
the largest eigenvalue of T+ — T~ , which determines the 
growth rate of the physical ground state of the system. 
The severity of the sign problem depends upon the dif- 
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ference between the largest eigenvalue of T+ + T~ and 
the largest eigenvalue of T+ — T~ — and thus on the 
prevalence of negative off-diagonal elements of T — and 
upon the concentration of psips required to achieve a suf- 
ficient rate of annihilation for the growth of the in-phase 
signal to be suppressed. A difficult sign problem does not 
necessarily imply a strongly correlated system. 

Eq. [2T] models the time dependence of the populations 
of positive and negative psips in a simulation in which all 
annihilation is forbidden. In a real FCIQMC simulation, 
however, the annihilation rate per psip does not tend to 
zero as the psip density tends to zero. The probability 
of encountering and annihilating an unrelated psip van- 
ishes, but a parent may still spawn a child of the opposite 
sign on its own determinant, in which case the parent and 
child annihilate each other. The low-density limit of the 
standard FCIQMC algorithm is better modeled by defin- 
ing T+ and T~ in a different way: T+ now contains all 
diagonal elements of T, regardless of sign, and all pos- 
itive off-diagonal matrix elements; whilst T~ is zero on 
the diagonal but contains the absolute values of the nega- 
tive off-diagonal elements. This redefinition corresponds 
to a change in viewpoint: instead of allowing parents to 
spawn children of the opposite sign on their own determi- 
nant, with subsequent annihilation, the negative diagonal 
elements of T+ remove psips from the simulation in one 
step, introducing an exponential decay of the psip pop- 
ulation on determinants for which < 0. The psip 
densities n'^ and nr remain positive at all times and 
the above analysis is unchanged, but the severity of the 
sign problem, as measured by the difference between the 
largest eigenvalues of T+ + T~ and T+ — T~ , is smaller 
than the above analysis suggests. 

From now on it will be assumed that T+ and T~ are 
defined as explained in the previous paragraph, and hence 
that all diagonal elements of T~ are zero. 



IV. POPULATION DYNAMICS 



where a = cq - n(0). Hence the evolution equation for p 
becomes 

^ = J2 ^iPi - «Pi + ««'e2^»-^c2i. (29) 



Eq. [29] is difficult to solve exactly, but the population 
dynamics it embodies can be illustrated using a simple 
one-component analogue: 



dp 



(30) 



where p{r) is the total psip population at time r, Fmax 
is larger than Tmax, and n(r) = noe^"^^^'^. It is common 
to start an FCIQMC simulation with a population of 
positive psips only, in which case the initial condition is 
p{0) = n(0) = no- As the initial psip population and an- 
nihilation rate are small, p(r) grows exponentially at the 
start of the simulation: p(r) = noe^"'^'''^. The exponen- 
tial growth ceases when p ^ Vmax/i^^ at which point the 
Hip'^ annihilation term (which is larger than the Hcn'^ term 
because Vmax > ^max) counteracts it. The psip popula- 
tion then enters a plateau phase, remaining rougly con- 
stant until the n(r) signal, which has been steadily grow- 
ing like noe-^"^^^'^, begins to dominate. From then on the 
population rises exponentially again, although now at a 
rate determined by Tmax- The ground-state wave func- 
tion thus spontaneously emerges from the simulation. 
The plateau begins to appear at a time ri determined 



by the equation p{ri) ^ uqc^ 



Kiax/^. Hence 



ri 



ln(14iax/^nQ) 

^^^ax 



(31) 



More precisely, solving Eq. [30] with the tzn'^ term omitted 
and assuming that i^no/Vmax ^ 1 (which must be the 
case if the annihilation rate is negligible at the beginning 
of the simulation) shows that p{r) reaches 95% of its 
plateau value at time 



Separating the psip population into positive and neg- 
ative contributions also explains the population dynam- 
ics observed during an FCIQMC simulation. Let V = 
T+ + T~, p = n+ + n~ and n = n+ — n~. The decou- 
pled ODEs in Eq. [26] can thus be written as 



dpi 



dr 

dn\ 
dr 



(27) 



In the long-time limit, n(r) becomes dominated by its 
ground-state component, which grows at a rate deter- 



^95% 



ln(19(Kiax/^no - 1)) 

Vmax 



(32) 



The end of the plateau occurs at a time T2 determined 



by the equation n(r2) = riQe^ 



"^max/^, and hence 



r2 



ln(14iax/^no) 



Combining Eqs. [31] and [33] yields 



(33) 



(34) 



mined by the largest eigenvalue, Tj^ 
transition matrix T: 



n(r) 



ae 



Co, 



The energy shift S appears on the diagonal of T = 
S — Eq^ of the S*! — H and is incorporated into the diagonal elements of 
T+; it therefore affects T = T+ - and V = T+ + 
equally. This dependence can be made explicit by writing 
(28) T^ax = ^+To and Kiax = ^ + where Tq and Fo(> Tq) 



8 



6.0 - 



S 4.0 

Ph 
O 



2.0 



0.0 



S-- 
S-- 
S-- 



Hoo + 0.5t 
ifoo - 0.5t 



0.0 



5.0 



10.0 

Ur 



15.0 



20.0 



FIG. 4. Effect of the initial value of the shift on the plateau 
height in an FCIQMC simulation of the 2D 18-site square 
lattice half- filled Hubbard model at U = 2t and k = (0,0). 
Periodic boundary conditions are applied to a 45° tilted cell 
with sides 3v^. The reference determinant, Do, is formed 
from occupying the 18 lowest-energy Bloch functions. Setting 
the shift slightly below Hqq = {Do\H\Do) reduces the plateau 
height but increases the time spent in the plateau region. 



are the largest eigenvalues of T and V when S 
leads to the equation 

n ^ S^To ~ S-Eo' 



0. This 



(35) 



In most FCIQMC simulations the shift is initially held 
at a fixed value, typically the Hartree-Fock energy, until 
the desired psip population is reached. If S is reduced 
and allowed to approach Eq from above, the simulation 
never exits the plateau. The duration of the plateau can 
be shortened by moving the initial value of S further 
above the ground-state energy Eq, but only at the cost of 
increasing the psip population {S -\-Vo)/k, at the plateau. 
The effect of modifying the initial shift in an FCIQMC 
simulation is shown in Fig. HJ 

The ODE in Eq. [30]is an example of a Riccati differen- 
tial equation^4 ^j^^ ^an be transformed from a quadratic 
first-order ODE into a linear second-order ODE using the 
substitution 



. . 1 du 
tvU dr 

The solution of the resultant second-order ODE is 



(36) 



u{t) =ci • oi^i ; 1 



2T^ 



^^max/2T'rn 



■ oi^i ; 1 



2T^ 



(37) 



where 



A^V(0)e 



2T^ 



4T2 

max 



(38) 



ci and C2 are constants of integration, o^^i is a confluent 
hypergeometric limit function,^^ and we have assumed 
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FIG. 5. Model population dynamics in an FCIQMC simula- 
tion. The psip population evolves according to Eq. [30] and 
shows the key features of the population dynamics in real 
FCIQMC calculations: an initial exponential growth phase 
followed by a plateau followed by a second slower exponential 
growth phase. As expected, the psip population reaches 95% 
of its plateau value when VmaxT95% ~ ln(19(ymax/^i^^o — 1)), 
the height of the plateau is given by p/no = Vmax//^no, and 
the plateau ends when Kiaxr2 ln(ymax//^no)Vmax/Tmax. 
The vertical dashed (dotted) lines show VmaxT95% (VmaxT2) 
for each set of parameters. The hypergeometric function, oi^i, 
was calculated using Ref. ^. 



that l^ax/ (2Tmax) IS uot an integer. The normalization 
of u(t) drops out of p(r) = {i<iu)~^du/dT^ leaving one 
arbitrary constant to be fixed by the initial conditions. 
Fig. [5] shows the resulting population dynamics for three 
different sets of parameters, imposing the boundary con- 
dition p(0) = no in all cases. Whilst this is an extremely 
simple model, it captures all of the features of the popu- 
lation dynamics seen in actual FCIQMC calculations. 

Let us now return to Eq. [291 and consider what happens 
beyond the end of the plateau, where the psip population 
begins to rise again until the energy shift is adjusted and 
a steady state with dp\/dt = is attained. Since the lin- 
ear V term is negligible in comparison with the quadratic 
terms in this regime, Eq. [29] becomes 



= ^ 
dt 



(39) 



implying that pi ^ ±ni and hence that ^ or ^ 
0. Determinants on which the ground- state amplitude 
coi is positive are occupied only by positive psips and 
determinants on which coi is negative are occupied only 
by negative psips. The signed psip density rii = — 
is proportional to coi and the unsigned density Pi = n'^ -\- 
nr to |coi|. 

We can also use Eq. [29]to understand why the plateau 
psip populations plotted in Fig. [2] are proportional to the 
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Hubbard U. Assuming that the growth of the in-phase 
signal, p, is fast enough to allow the plateau to emerge 
before the ground-state signal becomes significant, the 
plateau occurs when 

Y.^m^'^T.Pi- (40) 

ij i 

If U/t is large, the kinetic energy contributions to the 
Hamiltonian matrix (and hence T and V) can be ne- 
glected and changing U simply scales these matrices. 
Defining Vij = UV^'^ and p\ = Up[ leads to the follow- 
ing condition for the emergence of the plateau: 

E^jpj-^E^'S'. (41) 

ij i 

where V and contain no dependence upon U. The 
total psip population at the plateau, XliPi — ^^iPi^ 
therefore scales linearly with U. 

V. SIGN-PROBLEM-FREE SYSTEMS 

If there exists a similarity transform that maps T into 
V, and hence makes every off-diagonal element of T pos- 
itive and every off-diagonal element of H = 51 — T neg- 
ative, then V and T have identical eigenvalues. Even 
without annihilation, the in-phase and out-of-phase sig- 
nals, p and n, grow at the same rate and there is no 
difficulty extracting the ground-state (out-of-phase) sig- 
nal n = n+ — n~ . Such systems are sign-problem free and 
FCIQMC simulations of them yield correct ground-state 
energies with arbitrarily small psip populations, although 
no plateau phase is seen. Increasing the psip population 
serves only to reduce the stochastic error. Indeed, by for- 
mulating the problem in the transformed basis, we can 
carry out an FCIQMC simulation in which no negative 
psips need ever appear. 

A particularly simple type of similarity transformation 
does no more than change the signs of some of the basis 
determinants. No choice of signs is sufficient to render 
all off-diagonal elements of Ty positive in most cases, but 
there are a few interesting model systems in which sim- 
ple sign-changing transformations are effective. In some 
models, for example, all psips spawned on any given de- 
terminant Di have the same sign regardless of the loca- 
tion of their parent, so that positive and negative psips 
never mix. A simple sign-changing transformation that 
multiplies every basis determinant by the sign of the psips 
that occupy it then makes all off-diagonal elements of T 
positive. 

The antiferromagnetic Heisenberg model defined by 
the Hamiltonian 

H = J^Sr -Srs (42) 

where J>0, Sr is the vector spin operator on lattice site 
r, and the sum runs over nearest neighbors, is such a sys- 
tem if the lattice is bipartite. (Note that the basis states 



in this example are spin configurations rather than Slater 
determinants.) Indeed, the sign structure of the ground- 
state wave function of the bipartite Heisenberg model has 
long been known*^ Every non-zero off-diagonal matrix 
element of the Heisenberg Hamiltonian flips a neighbor- 
ing pair of spins from down-up to up-down or vice- versa. 
The Hamiltonian matrix element is positive if J > 0, 
implying that the corresponding matrix element of T is 
negative, so the signs of children produced by off-diagonal 
spawning events always differ from the signs of their par- 
ents. Initially it appears as if this ought to cause a seri- 
ous sign problem. Consider, however, a bipartite system 
consisting of two inter-penetrating sub-lattices, arranged 
such that the pairs of spins flipped by off-diagonal ele- 
ments of H are always on different sub-lattices. The ac- 
tion of any off-diagonal element of H increases the total 
value of Sz on one sub-lattice by 1 and decreases the total 
value of Sz on the other sub-lattice by 1. This allows all 
spin configurations to be assigned to one of two classes, 
A or B. Starting from the classical Neel state, which is 
arbitrarily assigned to class A, all states that differ from 
the Neel state by an even number of applications of H 
are also assigned to class A, while states that differ by 
an odd number of applications of H are assigned to class 
B. If an FCIQMC simulation is initialized by placing a 
population of positive psips on the Neel state, all psips 
created on configurations in class A will be positive while 
all psips created on configurations in class B will be neg- 
ative. Psips of different signs will never mix, no cance- 
lation is required, and there is no sign problem. In fact, 
by applying a simple unitary transformation in which the 
sign of every configuration in class B is flipped, all off- 
diagonal matrix elements of H can be made negative*^ 

By way of an example, the Heisenberg model for a 
two-dimensional 6x6 square lattice with periodic bound- 
ary conditions has a Hilbert space of 9.08 x 10^ con- 
figurations. Using just 1.8 X 10^ psips we find the 
E/N = —0.67886(2) J, which is in excellent agreement 
with the value E/N = -0.678871(8) J obtained by Runge 
using Green's function Monte Carlo. Note that, unlike 
Runge, we did not use importance sampling. FCIQMC 
simulations of Heisenberg models defined on other (non- 
bipartite or "frustrated") lattices display the same be- 
havior as for fermionic systems (Fig. [6]); in particular, 
the population plateau is a universal feature unless the 
system is sign-problem free. 

The Hubbard Hamiltonian, Eq. [191 is closely related 
to the Heisenberg Hamiltonian and is also sign-problem- 
free, although only in one dimension. Consider a ID 
Hubbard lattice with Ng sites and (N^) spin- up (spin- 
down) electrons, where the system is not necessarily half- 
filled and periodic boundary conditions are applied. If Ng 
is even and and are both odd or Ng is odd and 
and are both even, then there exists an analagous 
transformation to that described above for the Heisen- 
berg model. (This is most easily seen if the orbitals are 
ordered first by spin and then by their position in the 
lattice.) The ground- state energies of large ID Hubbard 
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FIG. 6. The population dynamics of an FCIQMC simulation 
of the antiferromagnetic Heisenberg model for a 5 x 5 trian- 
gular lattice. No periodic boundary conditions were imposed. 
As with fermionic systems, frustrated Heisenberg models dis- 
play a population plateau. 



models with many different choices of and can 
therefore be found using the FCIQMC method with very 
small numbers of psips. 

This brings us on to an important point: the sign prob- 
lem in FCIQMC is not constant for a given system but 
is dependent upon the choice of basis. Consider two 
Hamiltonian matrices, Hi and H2, which describe the 
same system but have been constructed using different 
basis sets, one of which is obtained from the other by a 
transformation matrix S, so that Hi = S~^H2S. The 
corresponding transition matrices Ti and T2 are related 
in the same way. However, when we partition the two 
transition matrices into T+ and T~, the partitions are 
not necessarily related by the same unitary transforma- 
tion. In other words, although Ti = 8^^X28, and hence 
- = S-^(TJ - T2")S, it is not in general the case 
that = S~^T^S. Consequently, the plateau height 
(and thus the ease with which FCIQMC can be used to 
find the ground state) depends on the basis in which the 
Hamiltonian matrix is expressed. 

An FCIQMC simulation of the 18-site half-filled ID 
Hubbard model at /7 = t exhibits no plateau when the 
Hamiltonian matrix is expressed in a basis of determi- 
nants of local orbitals (as in Eq. [19]), since this partic- 
ular problem is sign-problem free as explained above. 
Yet an FCIQMC simulation of exactly the same Hamil- 
tonian expressed in a basis of determinants of Bloch 
functions (as in Eq. [20]) exhibits a plateau at 6.9 x 10^ 
psips. As expected, the sign problem increases the dif- 
ficulty of the FCIQMC simulation in the Bloch repre- 
sentation. A short simulation using just 2.8 x 10^ psips 
in the local orbital basis gave a ground-state energy of 
— 18.8423(3)t, whereas 2.3 x 10^ psips were required to ob- 
tain a ground-state energy of — 18.84248(8)t in the Bloch 
orbital basis. Furthermore, our implementation currently 
conserves crystal momentum when using Bloch orbitals 
but does not make use of symmetry when using local 
orbitals; the Hilbert space used in the calculation with 
local orbitals therefore contains 2.36 x 10^ determinants. 



whereas that with Bloch orbitals contains only 1.31 x 10' 
determinants. 



VI. DISCUSSION 

The above observations provide some degree of insight 
into the FCIQMC method and the factors that determine 
how hard it is to apply FCIQMC to any given physical 
system. We now briefly comment upon other topics re- 
lated to FCIQMC before summarizing our work. 

The initiator approximation to FCIQMC (i-FCIQMC), 
whereby spawning events onto previously unoccupied de- 
terminants are only allowed if the parent determinant has 
a psip population above a specified threshold, has been 
shown to dramatically reduce the number of psips re- 
quired to obtain excellent results in many systems. — In 
i-FCIQMC, the effective Hilbert space grows and changes 
dynamically as the simulation proceeds, including only 
determinants that have a significant psip population or 
lie close (in terms of applications of the Hamiltonian) to 
determinants with a significant psip population. As a re- 
sult, psips are prevented from spawning in regions of low 
psip density and the annihilation rate is greatly increased 
for a given psip population. No plateau is observed in 
i-FCIQMC calculationsi^ii^ because the growth of the 
in-phase combination is suppressed by annihilation even 
when the total psip population is low. The i-FCIQMC 
approximation becomes increasingly good as the number 
of psips is increased because more of the Hilbert space 
becomes accessible and the psip dynamics more closely 
resembles that of the true Hamiltonian. 

In coupled cluster Monte Carlo^ (CCMC), the excita- 
tion amplitudes of the coupled cluster wave function are 
stochastically sampled in manner analogous to the sam- 
pling of configuration amplitudes in FCIQMC. Although 
the CCMC algorithm is much more complicated than the 
FCIQMC algorithm, we believe that a similar analysis 
holds. The particles or "excips"^^ that sample the dis- 
crete excitor space in CCMC also have positive or nega- 
tive signs and annihilation is crucial for the desired (and 
physically meaningful) solution to emerge. It is observed 
that CCMC calculations have a higher plateau than con- 
figuration interaction quantum Monte Carlo (CIQMC) 
calculations at the same truncation level^^ It is likely 
that this is because the excips have to explore a larger 
effective Hilbert space than the psips in the equivalent 
CIQMC simulation, so more excips are required to make 
the annihilation rate high enough to suppress the in- 
phase signal. 

In summary, we have shown that the sign problem in 
FCIQMC simulations is caused by the growth of an un- 
physical dominant solution of the coupled ODE's that 
describe the time evolution of the densities of positive 
and negative psips. This unphysical solution grows faster 
than the ground-state solution we seek and masks the 
ground-state signal. A similar problem arises in DMC 
simulations of fermionic systems, where the dominant 
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solution is the boson ground state. FCIQMC differs 
from DMC, however, in that the unphysical solution in 
FCIQMC is not an eigenfunction of the Hamiltonian op- 
erator. 

The annihilation of psips of opposite charge suppresses 
the growth of the unphysical dominant solution and 
allows FCIQMC simulations to converge on the true 
ground state, but works only when a minimum (and 
system-dependent) threshold in the psip population is 
exceeded. The combination of spawning and annihila- 
tion also leads to the characteristic population dynamics 
observed in FCIQMC calculations. The annihilation of 
psips may not be the most efficient way of suppressing the 
growth of the unphysical solution, and there remains a 
tantahsing possibility that an alternative approach could 
lead to a superior algorithm applicable to much larger 
systems. 
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